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Abstract 

We use asymptotic methods to study the early time stability of isotropic and homogeneous 
solutions filled with radiation which are close initially to the exact, flat, radiation solution 
in quadratic lagrangian theories of gravity. For such models, we analyze all possible modes 
, of approach to the initial singularity and prove the essential uniqueness and stability of the 

resulting asymptotic scheme in all cases except perhaps that of the conformally invariant Bach- 
Weyl gravity. We also provide a formal series representation valid near the initial singularity of 
the general solution of these models and show that this is dominated at early times by a form in 
which both curvature and radiation play a subdominant role. We also discuss the implications 
of these results for the generic initial state of the theory. 

^ ■ 

X: 

1 Introduction 

Cosmology in higher order gravity theories has become a field of special importance for a num- 
ber of basic questions about the structure of the early universe. In particular, stability issues 
for cosmological solutions in higher order gravity theories currently form a subfield of general- 
ized cosmological theory, one that combines interesting mathematical problems with the search for 
fundamental physical effects that might have taken place during the earliest moments of the uni- 
verse. The cosmological stability problem has two main facets, the perturbation theory aspect used 



* email: skot@aegean.gr 

' email: gkolionis@aegean.gr 

' email: atsok@aegean.gr 



1 



mainly in studies of structure formation, and the asymptotic stability aspect, basic for nonlinear 
and geometric dynamics issues. The latter problem naturally emerges when one is interested to see 
whether or not a given exact solution of the field equations has possibly a more general significance, 
whether a certain set of solutions have particular properties, or whether one needs to decide what 
can happen to the fate of any particular universe in this context if one waits long enough. 

The problem of asymptotic stability of isotropic and homogeneous cosmological solutions in 
higher order gravity theories has itself two major aspects, late-time stability-thai is deciding the fate 
of these universes in the distant future, and early-time stability-examining the past evolution at early 
times towards a possible initial singularity. Renewed interest in late time evolution cosmological 
problems in higher order gravity emerged after the seminal work of Barrow and Ottewill [1] where 
the issue of existence and stability of various cosmological solutions with an emphasis on the de 
Sitter and FRW ones was studied (see also [2] for related general late-time stability results for 
Friedman universes in the same context). Two important late-time stability problems that have 
received a lot of attention are the cosmic no-hair conjecture and the recollapse problem. For the 
original cosmic no-hair theorems in higher order gravity, see [3j[H[5], f° r limitations of this property 
cf. [U E] and for generalizations in higher order gravity of the corresponding situation in general 
relativistic cosmology cf. [U \10\ fTTj PT2l PT3l I14j . For various recollapse theorems in generalized 
cosmological theory see \15\ I16j. and for more elaborate and complete approaches and results, cf. 
|17l 1181 [T9] . It is interesting that the two issues, cosmic no-hair and recollapse of closed models are 
not unrelated, in fact a "premature" recollapse problem in closed universes that inflate has been 
formulated and studied, cf. [20J. 

The early-time evolution problem for isotropic and homogeneous cosmologies in higher order 
gravity is concerned with understanding the host of different behaviours that may occur near the 
initial singularity. The duality between bouncing and singular early time solutions was already 
present in the first papers of the subject, cf. |214 122j. These first solutions were shown to have the 
remarkable property of being horizon-breaking [23], but it was later pointed out that they may be 
unstable [24] . These first results also led to the realization that the question of the possible early 
time asymptotes of the admissible cosmological solutions of the higher order gravity equations was 
more involved [25], in fact, it was realized that even the simplest "radiation fluid" solution i 1 / 2 , 
being a solution in both the radiation filled case and in vacuum in these theories, may lead to a 
completely different set of properties than the situation encountered in general relativity. General 
properties of stability and conditions for instability of the early time flat and curved, radiation- 
filled, isotropic solutions were investigated in [26J, wherein although the interest then was in finding 
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instability properties of these systems as they are followed towards the initial singularity, there 
are various possible stability results one may deduce from the general conditions and equations 
by selecting appropriately various constants to take on specific forms. Stable solutions in vacuum 
near the initial singularity were obtained in the interesting works |28|, [29] when one includes in the 
basic quadratic lagrangian form R + aR 2 a term of the form Ric 2n ,ra £ Q, and considers a flat 
isotropic cosmology. The method of these papers was to use a linear perturbation analysis of the 
t 1 / 2 solution in vacuum and for flat FRW universes, and show that the various perturbations decay 
to zero asymptotically at early times. It is also known that this vacuum solution is stable under 
anisotropic, spatially homogeneous perturbations, cf. |30[ 13 lj . It is an interesting open problem 
to decide therefore the precise extent that a generic perturbation of the flat, vacuum, t 1//2 solution 
occupies in the whole space of solutions of the higher order gravity equations. 

Consequently, there are two separate early-time asymptotic problems concerning the stability 
of the flat t 1 / 2 FRW solution in higher order gravity: Its stability as a solution of the vacuum field 
equations, and that as solution of the higher order gravity field equations filled with a radiation 
fluid. Both problems need to be examined at all the various levels of stability: For the flat, radiation 
solution, we know [32] that it is stable in the space of all flat solutions of the theory asymptotically 
at early times in four spacetime dimensions. It remains to examine the precise behaviour of this 
solution consecutively with respect to curved FRW perturbations, anisotropic perturbations and 
generic inhomogeneous perturbations. As for the vacuum early time problem, as noted above, 
for the flat, vacuum solution, there are clear indications that it is stable with respect to various 
FRW and anisotropic perturbations, the strongest known results being its stability with respect to 
anisotropic perturbations [31], and with respect to perturbations in the R + (3R 2 action by adding 
a term of the form Ric 2 ™, cf. [281129]. 

In this paper, we focus on the problem of the early-time behaviour of the flat radiation t 1 / 2 

solution of higher order gravity with respect to curved FRW perturbations. That is, considered 

as a solution of the curved FRW equations for the R + (3R 2 action, what is the behaviour as we 

approach the initial singularity, i.e., as t — > of all solutions which are initially (that is, for some 

t* > 0) near this radiation solution? For this purpose, we approach the problem via the use of the 

method of asymptotic splittings developed in [331 134j . and trace all possible asymptotic behaviours 

that solutions to the higher order curved FRW equations may develop at early times. Following 

this geometric approach, we are able to show that the radiation solution is stable asymptotically 

at early times, meaning that the initial state of these universes proves to be a very simple one 

1 It is well known that the corresponding radiation solutions in general relativity are unstable with respect to any 
kind of perturbation and are also non generic, cf. [27] . 
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indeed. Given that this theory is known to admit an inflationary stage [35], this also means that 
any pre- inflationary period in such universes is necessarily isotropic and flat. 

The plan of this paper is as follows. In the next Section, we write down the field equations 
as a dynamical system in suitable variables thus defining a vector field, the curvature-radiation 
field. In Section 3, we find the possible dominant asymptotic forms which might be assumed by 
the curvature-radiation field towards the initial singularity, and prove the essential uniqueness of 
the acceptable asymptotic decompositions. In Section 4, we construct formal series expansions 
compatible with the dominant balances at the initial state and prove that these developments 
correspond to a local expansion of the general solution around the singularity. We conclude in 
Section 5, also commenting on our current results in the larger picture of ongoing and possible 
future work. 



2 The basic curvature-radiation field 

Below we shall be interested in the nature of the dynamics of FRW universes near the initial 
singularity. These cosmologies are determined by the Robertson- Walker metric of the form 

9i = -dt 2 + a 2 g 3 , (2.1) 

where a(t) denotes the scale factor, while each slice is given the 3-metric 

93 = V ^ r ^dr 2 + r 2 g 2 , ( 2 - 2 ) 
1 — kr A 

k being the (constant) curvature normalized to take the three values 0, +1 or —1 for the complete, 
simply connected, flat, closed or open space sections respectively, and the 2-dimensional sections 
are such that 

g 2 = d9 2 + sin 2 9 d(p 2 . (2.3) 

We assume that these spaces are filled with a radiation fluid with energy-momentum tensor T^ u = 
(p + p)UpU v + pg^ u , where the fluid velocity 4-vector is = 5q, and we take the equation of state 
to be p = p/S. 

Our general higher order action is I 2 ] 

S = - / AotaJtfyig, (2.4) 
2 Jm 4 

where Aotal is the lagrangian density of the general quadratic gravity theory given in the form 
Aotal = + ^matter j with 

C{R) = J R + /3i? 2 + 7Ric 2 + JRiem 2 , (2.5) 



"We set 8nG = c = 1, and the sign conventions are those of [36] 
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where /3,7, 8 are constants. Since in four dimensions we have the Gauss-Bonnet identity, 

5 R 2 GB d^ g = 0, R 2 GB = R 2 -4Ric 2 + Riem 2 , (2.6) 

Jm 4 

in the derivation of the field equations through a g- variation of the action associated with (|2.5p . 
only terms up to Ric 2 will matter. Therefore the variational derivative of the action leads to the 
following field equations: 



2RR> AU - -R 2 g pv - 2{g^g vu - g pu g pa )V p V a R 



+ 7 



2R pp R v p - -g pu R pa R pa + V P V p R pu + g pu V p V a R ap - 2V P V V R PP 



(2.7) 



For a homogeneous and isotropic space, we have a second useful identity, 

5 I (R 2 - 3Ric 2 )d^ s = , (2.8) 
Jm 4 

which further enables us to include the contribution of the Ric 2 term into the coefficient of R 2 , 
altering only the arbitrary constants. Hence, the field equations (12.7|) simplify as follows: 



2 y 6 



2RR pu - l -R 2 g pv - 2(g» p g va - g^g pa )V p V a R 



rpfJLV 



(2.9) 



where we have set £ = 2(3/3 + 7). This naturally splits into 00- and ii-components (i = 1, 2, 3), but 
only the 00-component of (12. 9p will be used below. 

Using the metric \2.1\ the field equation (|2.9[) leads to our basic cosmological equation in the 
form 



k + d l 



a a da 2 d 2 a . 
2 =- + ^ 5 2 ~~ 3 —t — 2k— 7 -| j 



d 2 k 21 



where £ is a constant defined by the constraint 

p e 



, (from V /t T^° = 0). 



(2.10) 



(2.11) 



3 a 4 

Setting x = a, y = a and z = d, Eq. (|2.10j) can be written as an autonomous dynamical system of 
the form 

x = ffc,RAD(x), x={v,y,z), (2.12) 
that is we have the dynamical system 



x 

y 

z 



y, 



(2.13) 



c 



fc 2 £ 3y 3 z 2 yz y k ky 
2£x 2 y + 2^2 + 2^~~~2£~2^ + ^2' 
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equivalent to the curvature-radiation vector field f & rad : M 3 — s- M 3 : (x,y,z) i— > fk,RAD( x > Vi z ) 
with 

/ £ 2 — 3y 3 z 2 yz y k ky\ , .. 

f fc.RAD (a?, !/,*) = a, y, -^i — + ^ + o JJt ~ «~ + ~2 J • ( 2 - 14 

\ 2£x z y 2x z 2y x 2£ 2£y x z J 

The curvature-radiation field ffc,RAD> or equivalently the dynamical system ([2.13p : combines the 
effects of curvature and radiation and describes completely the dynamical evolution of any radiation- 
filled FRW universe in higher order gravity. In the following Section we shall see how this field 
can split asymptotically and determine all dominant modes that it can develop on approach to the 
initial singularity in higher order gravity. 



3 Asymptotic splittings of the curvature-radiation field 

We are particularly interested below in the behaviour of the universe described by ([2.140 near the 
initial singularity, taken at t = Qj. Any such initial state will be completely described by giving 
the possible modes of approach of the various solutions to it. These modes are in turn determined 
by the behavior of ([2.140 near the initial singularity. How can this behaviour be found? For this 
purpose, we shall use below the method of asymptotic splittings, cf. [331 134j . According to this 
method, the vector field ffc RAD is asymptotically decomposed in such a way as to reveal its most 
important dominant features on approach to the singularity. This leads to a detailed construction of 
all possible local asymptotic solutions valid in the neighborhood of the finite-time singularity. These 
provide in turn a most accurate picture of all possible dominant features that the field possesses 
as it is driven to a blow up (for previous applications of this asymptotic technique to cosmological 
singularities, apart from [32], we refer to [37J). 

We expect that the vector field ffc,RAD will show some dominant features as we approach the 
singularity at t = 0, and these will correspond to the different, inequivalent ways that it splits in 
the neighborhood of the blow up. To describe the situation precisely, we need two definitions. We 
say that a solution b(t) of the system ([2.130 is asymptotic to another solution a(t) provided that 
the following two conditions hold (the first is subdivided): 

(i) Either (1) a(t) is an exact solution of the system, or (2) a(t) is a solution of the system 
(substitution gives = 0) as t — > oo, 

(ii) b(t) = a(t)[l + g(t)], g(t) -t 0, as t -> oo. 



3 The position of the initial singularity is really arbitrary, and we could have placed it at any to and used the 
variable r = t — to instead of t. 
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If either of these two conditions is not satisfied, then b(t) cannot be asymptotic to a(t). Secondly, a 
solution of the system fl2. 13f) is called dominant near the singularity if, for constants a = (9, rj, p) € 
C 3 , and p = (p, q, r) £ Q 3 , it is asymptotic to the form 

x (r) = at p = {et p ,r}T q ,pf). (3.1) 

For any given dominant solution of the system ()2. 13j) near the singularity, we call the pair (a, p) a 
dominant balance of the vector field ffc,RAD- 

Dominant balances, and the corresponding asymptotic integral curves, characterize vector fields 
near their blow up singularities^. The vector field ffc,RAD itself is decomposed asymptotically into 
a dominant part and another, subdominant part: 

f _ f (°) I f( sub ) /„ 9 n 

tfc.RAD - Ifc.RAD + fc,RAD ' K 6 - 1 ) 

and such a candidate asymptotic splitting (or decomposition) needs to be checked for consistency 
in various different ways before it is to be admitted as such. By direct substitution of the dominant 
balance forms in our system (|2.13p , we look for the possible scale invariant solutions of the syste nfl 

* = 4°rad(x)- (3-3) 

How many admissible asymptotic decompositions does the vector field f^RAD possess on approach 
to the initial state at t = 0? We recall that when k = 0, that is we have a flat, radiation-filled FRW 
model, the vector field f o,rad has two admissible asymptotic solutions near the initial singularity, 
as shown in [32J: In the first family, all flat, radiation solutions are dominated (or attracted) at 
early times by the form a(t) ~ i 1 / 2 , thus proving the stability of this solution in the flat case. There 
is a second possible asymptotic form near the singularity in the flat case, a(t) ~ t, but this contains 
only two arbitrary constants and hence it corresponds to a particular solution of the theory (cf. 
[32]). 

When k ^ 0, and we have the present situation of a radiation-filled, curved family of FRW 
universes to follow asymptotically near the past singularity, the field ffc,RAD has more terms-those 
that contain k in (I2.14p -than in the flat case. Since we, as noted above, already have a precise 
picture of the asymptotic forms of the flat case, we can now study in isolation the combined effects 
of curvature and radiation alone asymptotically. A simple combinatorial calculation shows that 



4 by a solution with a finite-time singularity we mean one where there is a time at which at least one of its 

components diverges. We note that the usual dynamical systems analysis through linearization etc is not relevant 

here, for in that one deals with equilibria, not singularities. 

5 A vector field f is called scale invariant if f(ar p ) = r p ~ 1 f(a), for a more detailed treatment, cf. [33] . 
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ffc,RAD (or the basic system (|2.13|) ) can decompose precisely in 

= 255 (3.4) 















H 






: 






k 2 1 







different ways. Each one of these 255 different modes leads to an asymptotic splitting of the form 
(|3.2p . each one of which may contain many possible dominant balances and needs to be checked for 
admissibility. Any candidate asymptotic splitting of the form (|3.2p will be acceptable in principle 
provided the candidate subdominant part tends to zero asymptotically, that is it indeed behaves 
as a subdominant contribution to the dominant asymptotic form the field splits into (cf. (|3.2p ). 

It would be exceedingly cumbersome to list here the detailed results we have for each particular 
decomposition, so we prefer to describe them in a qualitative way instead. The first qualitative 
conclusion from the asymptotic decomposition analysis we performed is that the first term in the 
field ffc,RAD, 

z '^e^J' (3 - 5) 

cannot be further splitted in its two parts: Each one of the resulting terms in (|3.5p . if taken alone 
asymptotically, leads to 64 decompositions for which the subdominant part either grows or tends 
to a constant, rather than decays to zero as t — ¥ 0. Of the remaining 127 decompositions, 64 
contain the term (|3.5p while 63 do not contain it in their dominant fi°j^rj P ar t. Also, 15 of the 



127 decompositions of the field do not contain in their dominant part neither the term (|3.5p nor 
any of the terms containing the curvature k. It is interesting that of these 15 decompositions, 
7 do not contain the linear term — y/2£ in the dominant part, and among these there is only 
one that eventually leads to a fully acceptable (cf. below) dominant balance: namely, this is 
the decomposition that does not contain (in its dominant part) any of the terms containing the 
parameters £,k as well as the linear term —y/2£, but does contain all the other remaining three 
terms of the vector field (|2.14p . We therefore arrive at the interesting qualitative conclusion that in 
the higher order gravity asymptotics towards the initial singularity of the curved, radiation-filled 
FRW models, neither radiation not curvature play a dominant role! 

The rest 126 decompositions fail to lead to an acceptable dominant balance for one of the three 
different reasons: Either 

1. the dominant exponents in a balance (cf. (|3.ip ) cannot be determined as a solution of the 
dominant system, or 

2. the condition for the existence of a acceptable subdominant part 

a 'T-i = °' (3 - 6) 



is not satisfied, or 

3. they lead to a conflict with the arbitrariness of the position of the singularity (this translates 
into an invalid condition on the eigenvalues of the so-called Kovalevskaya matrix, cf. [33] and 
below). 

We conclude that the only acceptable asymptotic splitting of the vector field f ^°rad ^ s °f * ne f° rm 

f , _ f (0) , n (sub) 

Ifc,RAD — RAD r fc,RAD' 



fLU(x)=( 2 /,^^ + ^-^), (3.7) 



with with dominant part 

f(o) ^ - (.. z V + f! _ K£ 

and subdominant part 

•0 ub ) ^ _ / n n C 2 ~ _ J/_ fey 

2£x 2 y 2£ 2£y a; 



f&M = 0,0,^^-^-^ + 3). (3. 



A final comment about the asymptotic splittings of the field equations of the present Section is in 
order. It is interesting that the dominant part of the vector field given by Eq. (|3.7p is precisely the 
same as that of the field in the flat, radiation dominated case treated in [32] (see Eq. (16) in that 
paper). Their difference lies in the subdominant parts of the two cases, the curved one treated here 
and the flat case in [32]: Here the subdominant part given by ()3.8p contains precisely the terms of 
the vector field f orad (radiation and linear terms), plus the three curvature terms (those with a k 
in Eq. §Z7\M . 

We have completed in this Section the first part of our asymptotic analysis via the method 
of asymptotic splittings, that is we found all possible asymptotic systems on approach to the 
initial singularity. This process amounts to dropping all terms that are small, and replace exact 
by asymptotic relations (by this we mean using (13. 3D in conjunction with (|3,2p and (13 .6p instead 
of (|2.12p and (I2.14p ). This first part of the application of the method of asymptotic splittings 
allows us to conclude that there is essentially a unique such system, and we were able to extract 
some qualitative results about the behaviour of our basic vector field without actually solving the 
systems. In the next Section, we shall proceed to study the solutions of our asymptotic system 
through the processes of balance, subdominance and consistency. 



4 Stability of the asymptotic solutions 

We now look for the possible asymptotic solutions, asymptotic forms of integral curves of the 
curvature-radiation field ffcRADj that is we search for the dominant balances determined by the 
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dominant part f^RAD §i ven by Eq. (|3.7p . For this purpose, we substitute in the dominant system 
(x, y, z)(t) = f/uRAD the forms (|3.ip and solve the resulting nonlinear algebraic system to determine 
the dominant balance (a, p) as an exact, scale invariant solution. This leads to the unique curvature- 
radiation balance Bk,RAD € C 3 x Q 3 , with 

Bl , MD = (a,p)=((4-f),(i,4,-?)), (4.1, 

where 9 is a real, arbitrary constant. In particular, this means that the vector field f J^ad ^ s a 
scale invariant system, cf . [33|, [M] . 

Further, we need to show that the higher order terms (|3.8p in the basic decomposition of the 
curvature-radiation field (|3.2p are themselves weight-homogeneous with respect to the curvature- 
radiation balance (j4. 1 j) for this to be an acceptable one. To prove this, we first split the subdominant 
part (|3.8p by writing 

fffi D (x) = fguDto + fgLnW + fl 3 L D (x), (4-2) 

where 

f£UW= (0,0,§), fg UD (x)= (0,0,^^-|), fgL D (x)= (0,0,-A). (4.3) 

We can then find the required condition that it is indeed subdominant. Using the balance Bk,RAD 
defined by Eq. we find that (note that f" 1 = i~ 5 / 2 ) 



fc,RAD 



£ P -1 fc,RAD* 

^^^^ - I&adM* (4-5) 



4 2 Ld(" p ) 

j.p-1 ~~ 1 fc,RAD 



%^ = «&*,««•■ (4-6) 

Hence, taking the limit as t — > 0, we see that these fractions go to zero asymptotically provided 
that the forms fj^ RAD (a),i = 1,2,3, are all different from zero. This happens only when £ ^ 0, 
that is for all cases except when 3/3 + 7 = 0. We conclude that this result is true in all higher order 
gravity theories except perhaps the so-called conformally invariant Bach-Weyl gravit>o, cf. [38J. 
Since the subdominant exponents 

q (°) = < g« = 1 < = 2 < = 3, (4.7) 
are ordered, we conclude that the subdominant part (|3.8p is weight-homogeneous as promised. 



6 Apparently, this case needs a separate treatment altogether. 
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Let us now move on to the last part of our asymptotic investigations. For our asymptotic 
solutions to be valid, we further need to check their consistency with the overall approximation 
scheme we are using. This means that we need to construct a series representation of the asymptotic 
solutions valid locally around the initial singularity, so that it is dominated by the dominant balance 
solutions we have built so far. The degree of generality of these final series solutions will of course 
depend on the number of arbitrary constants in them. We expect such constants to appear in 
certain places in the final formal developments we are after and, as an intermediate step in the 
construction of the series expansions, we now calculate the precise positions where such constants 
would appear. 

As explained in [33], the arbitrary constants of any (particular or general) solution first appear 
in those terms in the asymptotic series solution whose coefficients Cj have indices i = gs, where g 
is a non-negative /C-exponent and s denotes the least common multiple of the denominators of the 
set of all subdominant exponents (|4.7p and those of all the /C-exponents with positive real parts (in 
our case, s = 2). These exponents are numbers belonging to the spectrum of an important matrix 
in asymptotic analysis, the so-called Kovalevskaya matrix given by 

K = H°RAD( a )-diag( P ). (4.8) 

Hence, the /C-exponents depend of the dominant part of the vector field as well as the dominant 
balance. In our case, the Kovalevskaya matrix is 



/ -1/2 1 \ 



K 



fc,RAD 



(4.9) 



1/2 1 
\ -1/2 5/4 1/2 / 

with spectrum 

spec(/C fc)RAD ) = {-1,0,3/2}, (4.10) 

and corresponding eigenvectors 

{(4, -2,3), (4,2,-1), (1,2,2)}. (4.11) 

The number of non-negative /C-exponents equals the number of arbitrary constants that appear in 
the series expansions. There is always the —1 exponent that corresponds to an arbitrary constant, 
the position of the singularity, and because the spec(/Cfc,RAD) in our case possesses two non-negative 
eigenvalues, the balance rad indeed corresponds to the dominant behaviour of a general solution 
having the form of a formal series and valid locally around the initial singularity. To find it, we 
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substitute the Puiseux series expansions 

oo oo oo 

x(t)=J2c L {tT+*, y{t)=J2 C ^'K Z(t) =J^C3i*5 - ', (4.12) 
i=0 i=0 i=0 

where cio = #,C2o = #/2,C3o = —0/A, in the original system (|2.13p and we are led to various 
recursion relations that determine the unknowns cu,C2i,c^i term by term. Further algebraic ma- 
nipulations lead to the final series representation of the solution in the form: 



x{t) _ , ^ - - t>» + c a f + - ^ ) + ■ ■ ■ . (4.13) 

The corresponding series expansions for y(t) and z(t) are given by the first and second time deriva- 
tives of the above expression respectively. 

As a final test for admission of this solution, we use Fredholm's alternative to be satisfied by 
any admissible solution. This leads to a compatibility condition for the positive eigenvalue 3/2 and 
the associated eigenvector: This condition has the form 

fC - iA Cj = 0, (4.14) 



s 

where I denotes the identity matrix, and we have to satisfy this at the j = 3 level. This gives the 
following orthogonality constraint, 

/ -2ci 3 + c 2 3 \ 



(1,2,2) 



-C23 + C33 
V -3C13 + |c 2 3 - C33 J 



0. (4.15) 



It is immediate to check that this is indeed satisfied, thus leading to the conclusion that (I4.13P 
corresponds indeed to a valid asymptotic solution around the singularity. We also see that setting 
k = in the series expansion we are led exactly to the form found for the flat, radiation case, cf. 
Eq. (21)0 of [32] . 

Our series (|4.13p has three arbitrary constants, 9, C13 and a third one corresponding to the 
arbitrary position of the singularity (taken here to be zero without loss of generality), and is 
therefore a local expansion of the general solution around the initial singularity. Since the leading 
order coefficients are real, by a theorem of Goriely and Hyde, cf. [33], we conclude that there is 
an open set of initial conditions for which the general solution blows up at the finite time (initial) 
singularity at t = 0. This proves the stability of our solutions in the neighborhood of the singularity. 



7 In that reference, 12£8 3 was mistakenly written as 24£# 3 



12 



5 Conclusions 



In this paper we have considered the behaviour of the flat, radiation-filled isotropic and homoge- 
neous cosmological models in the general quadratic gravity theory near the initial singularity. We 
have shown that asymptotically the effects of curvature and radiation on the global evolution be- 
come negligible and the initial state of these universes is basically described by the flat, vacuum i 1 / 2 
solution of these theories (with the possible exception of the solutions in the conformally invariant 
Bach-Weyl theory - this case needs to be analyzed separately). In particular, all curved radiation 
solutions tend to this solution asymptotically as we approach the initial state in these theories. 

There are precisely 255 ways that the basic curvature-radiation vector field of this problem 
can split and each one of these needs to be checked for admissibility for it could in principle 
contain various dominant asymptotic solutions. However, taking into account various asymptotic 
conditions that have to hold in order to have admissible solutions, we are left with only one possible 
asymptotic decomposition of the vector field near the initial state. Further, it turns out that this 
asymptotic form admits only one dominant balance, exact solution of the scale invariant system, 
and this is precisely the form that behaves as i 1//2 near the singularity. Using various asymptotic 
and geometric arguments, we were able to built a solution of the field equations in the form of a 
Puiseux formal series expansion compatible with all other constraints, dominated asymptotically 
by the t 1 ^ 2 solution and having the correct number of arbitrary constants that makes it a general 
solution. In this way, we can conclude that this exact solution is an attractor of all homogeneous 
and isotropic radiation solutions of the theory thus proving stability against such 'perturbations'. 

The impressive restrictions placed by the higher order field equations on the structure of the 
possible initial cosmological states of the theory imply that the initial state of these radiation 
universes probably resembles that of a vacuum, flat model since, as we showed here, the only 
possible mode of approach to the singularity is one in which both the curvature and radiation modes 
enter asymptotically only in the subdominant part of the vector field. In fact, an investigation of 
the asymptotic stability of the flat, vacuum initial state in these theories with respect curved FRW 
perturbations is currently under scrutiny, cf. [39| . If both states prove to be stable ones in a way 
along the lines of the present paper, whether or not the radiation or the vacuum initial states are 
in addition generic ones in the space of solutions of the higher order gravity equations is currently 
an open problem. 
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